# This code is part of LINKEQ.
#
# (C) Copyright LINKE 2023.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
#
# -*- coding: utf-8 -*-
# @Time     : 2024/11/12 16:54
# @Author   : deviludier @ LINKE
# @File     : LargeCircuitExecution.py
# @Software : PyCharm
'''
This file contains the class of large circuit execution, provides a method to execute large circuit by cutting them to smaller circuits
'''
import os
import sys


# import depency
import sympy as sp
import numpy as np
from quantumcircuit.circuit import QuantumCircuit as QC
from quantumcircuit.gate import *
from scipy.optimize import linprog
import copy
import qiskit
import qiskit_aer

from qiskit_aer import AerSimulator
from qiskit import QuantumCircuit, transpile
from qiskit.circuit import Parameter
from qiskit import QuantumRegister, ClassicalRegister
from qiskit.quantum_info import Pauli, SparsePauliOp, Operator, DensityMatrix
# from qiskit_aer.noise import depolarizing_error, phase_damping_error, amplitude_damping_error, pauli_error, thermal_relaxation_error
# from qiskit.primitives import Estimator
from backend.simulator import *
from qiskit.circuit import Gate as gt
sys.path.append("/home/normaluser/deviludier/LinkeQ_v0.2/backend")

dict = {'cnot':[('i','i',0.5),('i','x',0.5),('z','i',0.5),('z','x',-0.5)],
        'cz':[('rz-pi/2','rz-pi/2',0.5),('rz+pi/2','rz+pi/2',0.5),
              ('1','i',-0.5),('i','1',-0.5),
              ('1','rz-pi',0.5),('i','0',0.5),
              ('0','i',0.5),('rz-pi','1',0.5),
              ('0','rz-pi',-0.5),('rz-pi','0',-0.5)]}
max_len = 256
class Observable():
    def __init__(self, pauli_list, coef_list) -> None:
        self.observable = {}
        for i in range(len(pauli_list)):
            self.observable[pauli_list[i]] = coef_list[i]

def printQC(circ):
    for i in range(1,len(circ.gate_list)):
        print(circ.gate_list[i].name,circ.gate_list[i].get_qubits(), end='')
    print('\n')

def optimize_x(U , S):
    """
    Solve the optimization problem:
    Minimize sum(|x[k]|) subject to sum(x[k] * S[k]) = U.

    Parameters:
    U : ndarray
        Target matrix (2D array).
    S : list of ndarray
        List of matrices (2D arrays).

    Returns:
    x_opt : ndarray
        Optimal solution x.
    """
    U = np.column_stack((np.real(U),np.imag(U)))
    S = [np.column_stack((np.real(m),np.imag(m))) for m in S]

    n_matrices = len(S)  # Number of matrices in the list
    dim = U.shape  # Dimensions of U and each S[k]
    n_elements = dim[0] * dim[1]  # Number of elements in each matrix

    # Flatten U and S for linear programming
    U_flat = U.flatten()
    S_flat = [s.flatten() for s in S]

    # Linear programming formulation
    # Variables: x and z, where len(x) = len(z) = n_matrices
    c = [0] * n_matrices + [1] * n_matrices  # Objective: min sum(z[k])

    # Equality constraints: sum(x[k] * S[k]) = U
    A_eq = np.zeros((n_elements, 2 * n_matrices))
    for k in range(n_matrices):
        A_eq[:, k] = S_flat[k]
    b_eq = U_flat

    # Inequality constraints: z[k] >= x[k] and z[k] >= -x[k]
    A_ub = np.zeros((2 * n_matrices, 2 * n_matrices))
    for k in range(n_matrices):
        A_ub[k, k] = 1  # -x[k]
        A_ub[k, n_matrices + k] = -1  # +z[k]
        A_ub[n_matrices + k, k] = -1  # +x[k]
        A_ub[n_matrices + k, n_matrices + k] = -1  # +z[k]
    b_ub = [0] * (2 * n_matrices)

    # Solve the linear program
    result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds =(-2,2), method='highs')

    if not result.success:
        raise ValueError(f"Optimization failed: {result.message}")

    # Extract x from the result
    x_opt = result.x[:n_matrices]
    return x_opt

class largeCircuitExecution():
    def __init__(self, qc, observable, label) -> None:
        self.target_qc = qc
        self.observable = observable
        self.label = label
        self.circuit_list = None
        self.observable_list = None
        self.full_map = None
        self.coef_list = None
        self.observable_coef_list = None
        self.parameter_list = {}


    def qpd(self, circ):
        pass

    def partitionObservable(self):
        pass

    def partitionCircuits(self):
        '''
        1. 根据label寻找到需要切割的两比特门；
        2. 调用twoQubitGateDecomposition进行两比特门切割，得到返回的量子线路组合；
        3. 根据返回的量子线路组合，构造整体分割后的量子线路组合；
        '''
        pass

    def circuitExecution(self):
        pass


    def CirctoList(self):
        gate_list = self.target_qc.gate_list
        clist = []
        for i in range(1,len(gate_list)):
            name = 'i'
            if gate_list[i].name == 'X':
                clist.append([('x',gate_list[i].get_qubits())])
            if gate_list[i].name == 'Y':
                clist.append([('y',gate_list[i].get_qubits())])
            if gate_list[i].name == 'Z':
                clist.append([('z',gate_list[i].get_qubits())])
            if gate_list[i].name == 'H':
                clist.append([('h',gate_list[i].get_qubits())])
            if gate_list[i].name == 'CX':
                clist.append([('h',(gate_list[i].get_qubits())[-1:])])
                clist.append([('cz',gate_list[i].get_qubits())])
                clist.append([('h',(gate_list[i].get_qubits())[-1:])])
            if gate_list[i].name == 'CZ':
                clist.append([('cz',gate_list[i].get_qubits())])
            if gate_list[i].name == 'RX':
                clist.append([('rx',gate_list[i].get_qubits(), gate_list[i].get_para())])
                if gate_list[i].is_symbol:
                    p = Parameter(gate_list[i].get_para())
                    self.parameter_list[gate_list[i].get_para()] = p
            if gate_list[i].name == 'RY':
                clist.append([('ry',gate_list[i].get_qubits(), gate_list[i].get_para())])
                if gate_list[i].is_symbol:
                    p = Parameter(gate_list[i].get_para())
                    self.parameter_list[gate_list[i].get_para()] = p
            if gate_list[i].name == 'RZ':
                clist.append([('rz',gate_list[i].get_qubits(), gate_list[i].get_para())])
                if gate_list[i].is_symbol:
                    p = Parameter(gate_list[i].get_para())
                    self.parameter_list[gate_list[i].get_para()] = p
        return clist
    
    def CircuitCutting(self, str, parameter_value = []):
        circ = self.target_qc
        circList = self.CirctoList()
        n=circ.qubit_number
        # map origin qubits to subsystem A and B
        mapa = []
        mapb = []
        fullmap = []
        for i in range(len(str)):
            if str[i] == 'a' or str[i] == 'A':
                mapa.append(i)
                fullmap.append(len(mapa)-1)
            if str[i] == 'b' or str[i] == 'B':
                mapb.append(i)
                fullmap.append(len(mapb)-1)
        self.full_map = fullmap
        #print(mapa,mapb,fullmap)
        #finding gates needed to cut and then cut them
        res = []
        subrout = []
        ce = []
        clist = []
        for i in range(len(circList)):
            gates = circList[i]
            assert(len(gates)==1)
            act_qubit = gates[0][1]
            lista = []
            listb = []
            for j in range(len(act_qubit)):
                if str[act_qubit[j]] == 'a' or str[act_qubit[j]] == 'A':
                    lista.append(act_qubit[j])
                if str[act_qubit[j]] == 'b' or str[act_qubit[j]] == 'B':
                    listb.append(act_qubit[j])
            if len(lista) == 0 or len(listb) == 0:
                continue
            clist.append(i)
            part = dict.get(gates[0][0])
            if part == None:
                raise ValueError(f'No qpd for {gates[0][0]} gate')
            subgate = []
            coef = []
            for v in part:
                coef.append(v[2])
                subgate.append([(v[0],[act_qubit[0]]) , (v[1],[act_qubit[1]]) ] )
            subrout.append(subgate)
            ce.append(coef)
        num = 1
        for v in subrout:
            num*=len(v)
        subcirc = []
        subcoef = []
        for i in range(num):
            ii = i
            m = num
            c = 1.
            subcirc.append(copy.deepcopy(circList))
            for j in range(len(subrout)):
                m = int(m/len(subrout[j]))
                i0 = int(ii/m)
                #print(i0)
                c *= ce[j][i0]
                subcirc[i][clist[j]] = subrout[j][i0]
                #print(subrout[j][i0])
                ii = ii%m
            subcoef.append(c)
        #print(subcoef)
        #print(subcirc)
        #generate sub circuit list
        suba = []
        subb = []
        for v in subcirc:
            cira = QuantumCircuit(len(mapa)) 
            qubits_a = QuantumRegister(len(mapa))
            qubits_b = QuantumRegister(len(mapb))
            cl_a = ClassicalRegister(len(mapa))
            cl_b = ClassicalRegister(len(mapb))
            clm_a = ClassicalRegister(max_len)
            clm_b = ClassicalRegister(max_len)
            ca = qiskit.QuantumCircuit(qubits_a,cl_a,clm_a)
            cb = qiskit.QuantumCircuit(qubits_b,cl_a,clm_b)
            circs = [ca,cb]
            q = [qubits_a,qubits_b]
            clm = [clm_a,clm_b]
            tag = [0,0]
            #cirb = QuantumCircuit(len(mapb))
            for gates in v:
                for g in gates:
                    #gate = Gate(g[0])
                    flag = 1
                    if str[g[1][0]] == 'a' or str[g[1][0]] == 'A':
                        flag = 0
                    if g[0] == 'i':
                        continue
                    if g[0] == '0':
                        #op = pauli_error([["I", 0.5], ['Z', 0.5]])
                        #print(op)
                        assert(tag[flag] < max_len)
                        circs[flag].measure(q[flag][fullmap[g[1][0]]], clm[flag][tag[flag]])
                        tag[flag] += 1
                        continue
                    if g[0] == '1':
                        #op = pauli_error([["I", 0.5], ['-Z', 0.5]])
                        #print(op)
                        assert(tag[flag] < max_len)
                        circs[flag].x(fullmap[g[1][0]])
                        circs[flag].measure(q[flag][fullmap[g[1][0]]], clm[flag][tag[flag]])
                        circs[flag].x(fullmap[g[1][0]])
                        tag[flag] += 1
                        continue
                    if g[0] == 'rz-pi/2':
                        circs[flag].rz(-np.pi/2,fullmap[g[1][0]])
                        continue
                    if g[0] == 'rz+pi/2':
                        circs[flag].rz(np.pi/2,fullmap[g[1][0]])
                        continue
                    if g[0] == 'rx-pi/2':
                        circs[flag].rx(-np.pi/2,fullmap[g[1][0]])
                    if g[0] == 'rx+pi/2':
                        circs[flag].rx(np.pi/2,fullmap[g[1][0]])
                    if g[0] == 'rz-pi':
                        circs[flag].rz(-np.pi,fullmap[g[1][0]])
                    if g[0] == 'rx':
                        circs[flag].rx(g[2],fullmap[g[1][0]])
                    if g[0] == 'ry':
                        circs[flag].ry(g[2],fullmap[g[1][0]])
                    if g[0] == 'rz':
                        circs[flag].rz(g[2],fullmap[g[1][0]])
                    if g[0] == 'x':
                        #gate = X(fullmap[g[1][0]])
                        circs[flag].x(fullmap[g[1][0]])
                    if g[0] == 'y':
                        #gate = Y(fullmap[g[1][0]])
                        circs[flag].y(fullmap[g[1][0]])
                    if g[0] == 'z':
                        #gate = Z(fullmap[g[1][0]])
                        circs[flag].z(fullmap[g[1][0]])
                    if g[0] == 'h':
                        #gate = H(fullmap[g[1][0]])
                        circs[flag].h(fullmap[g[1][0]])
                    if g[0] == 'cnot':
                        #gate = CX(fullmap[g[1][0]],fullmap[g[1][1]])
                        circs[flag].cx(fullmap[g[1][0]],fullmap[g[1][1]])
                    if g[0] == 'cz':
                        #gate = CX(fullmap[g[1][0]],fullmap[g[1][1]])
                        circs[flag].cz(fullmap[g[1][0]],fullmap[g[1][1]])
                    #if str[g[1][0]] == 'a' or str[g[1][0]] == 'A':
                        #cira.add_gate(gate)
                    #if str[g[1][0]] == 'b' or str[g[1][0]] == 'B':
                        #cirb.add_gate(gate)
            suba.append(circs[0])
            subb.append(circs[1])
            #print(circs[0].num_qubits)
        self.circuit_list = [suba,subb]
        self.coef_list = subcoef
        return suba,subb,subcoef
    
    def ObservableCutting(self, str):
        mapa = []
        mapb = []
        fullmap = []
        for i in range(len(str)):
            if str[i] == 'a' or str[i] == 'A':
                mapa.append(i)
                fullmap.append(len(mapa)-1)
            if str[i] == 'b' or str[i] == 'B':
                mapb.append(i)
                fullmap.append(len(mapb)-1)
        #print(mapa,mapb,fullmap)
        observables = self.observable.observable
        pauli_a = []
        pauli_b = []
        coef = []
        for k,v in observables.items():
            #print(k,v)
            pauli_stra = list('I'*len(mapa))
            pauli_strb = list('I'*len(mapb))
            i = 0
            while i<len(k):
                assert(k[i]=='X' or k[i]=='Y' or k[i]=='Z' or k[i]=='x' or k[i]=='y' or k[i]=='z' or k[i]=='i' or k[i]=='I')
                temp = k[i]
                i += 1
                assert(k[i]==' ')
                i += 1
                index = 0
                while i<len(k) and k[i] >= '0' and k[i]<='9':
                    index = 10*index+int(k[i])
                    i += 1
                assert(i>=len(k) or k[i]==' ')
                i += 1
                temp += f' {fullmap[index]} '
                if str[index] == 'a' or str[index] == 'A':
                    pauli_stra[len(mapa)-1-fullmap[index]] = temp[0]
                if str[index] == 'b' or str[index] == 'B':
                    pauli_strb[len(mapb)-1-fullmap[index]] = temp[0]
            pauli_a.append(''.join(pauli_stra))
            pauli_b.append(''.join(pauli_strb))
            coef.append(v)
        # Oa = Observable(pauli_a,coef)
        # Ob = Observable(pauli_b,coef)
        #print(pauli_a)
        self.observable_list = [pauli_a,pauli_b]
        self.observable_coef_list = coef
        return 
    
    def Execute(self,shots=1024):
        suba = self.circuit_list[0]
        subb = self.circuit_list[1]
        Oa = self.observable_list[0]
        Ob = self.observable_list[1]
        coef = self.coef_list
        assert(len(suba)>0)
        if len(self.parameter_list)>0:
            return
        da = suba[0].num_qubits
        db = subb[0].num_qubits
        o_coef  = self.observable_coef_list
        #print(Oa,Ob,coef,o_coef)
        result = 0.
        for i in range(len(coef)):
            circ_a = suba[i]
            circ_b = subb[i]
            #print(circ_a)
            #print(circ_b)
            result_i = 0.
            for j in range(len(o_coef)):
                result_a = 1.
                result_b = 1.
                if all(bit == 'i' or bit == 'I' for bit in Oa[j]):
                    result_a = 1.
                    simulator = AerSimulator()
                    # Run and get counts
                    circ_a.measure(circ_a.qregs[0][0],  circ_a.cregs[0][0])
                    results = simulator.run(circ_a,shots=shots).result()
                    counts = results.get_counts(circ_a)
                    sum = 0
                    for k,v in counts.items():
                        cr1_part = k[:max_len]          # 前 m 位为 cr1
                        cr2_part = k[-1:]       # 后 n 位为 cr2
                        # 检查 cr1 是否全为0
                        if all(bit == '0' for bit in cr1_part):
                            sum += v
                    result_a = float(sum)/float(shots)
                else:
                    mat_rhoa = np.zeros((1<<circ_a.num_qubits, 1<<circ_a.num_qubits),dtype= complex)
                    mat_rhoa[0,0] = 1
                    circ_a_exe = circ_a
                    oa_str = list(Oa[j])
                    c2_num = 0
                    tag = 0
                    for o in range(len(oa_str)):
                        if oa_str[len(oa_str)-1-o] == 'x' or oa_str[len(oa_str)-1-o] == 'X':
                            circ_a_exe.h(o)    
                            oa_str[len(oa_str)-1-o] = 'Z'
                        if oa_str[len(oa_str)-1-o] == 'y' or oa_str[len(oa_str)-1-o] == 'Y':
                            circ_a_exe.s(o)
                            circ_a_exe.h(o)
                            circ_a_exe.x(o)
                            oa_str[len(oa_str)-1-o] = 'Z'
                        if oa_str[len(oa_str)-1-o] == 'i' or oa_str[len(oa_str)-1-o] == 'I':
                            continue
                        circ_a_exe.measure(circ_a_exe.qregs[0][o],  circ_a_exe.cregs[0][tag])
                        tag += 1
                        c2_num += 1
                    #rhoa = qiskit.quantum_info.DensityMatrix(mat_rhoa)
                    #rhoa = rhoa.evolve(circ_a_exe)
                    # print(oa_str,circ_a_exe)
                    #opa = SparsePauliOp.from_list([(''.join(oa_str), flag)])
                    #print(rhoa)
                    #result_a = rhoa.expectation_value(opa).real
                    #print(rhoa)
                    #estimator = Estimator()
                    #result_a = estimator.run(circ_a, opa).result().values.real[0]
                    simulator = AerSimulator()
                    #print(circ_a_exe)
                    # Run and get counts
                    results = simulator.run(circ_a_exe,shots=shots).result()
                    counts = results.get_counts(circ_a_exe)
                    c1_num = 256
                    total_sum = 0
                    sum = 0
                    cr2_counts = {}
                    #print(counts)
                    for k,v in counts.items():
                        cr1_part = k[:c1_num]          # 前 m 位为 cr1
                        cr2_part = k[-c2_num:]       # 后 n 位为 cr2
                        # 检查 cr1 是否全为0
                        if all(bit == '0' for bit in cr1_part):
                            num_ones = cr2_part.count('1')
                            cr2_counts[cr2_part] = cr2_counts.get(cr2_part, 0) + v
                            parity = 1 if (num_ones % 2 == 0) else -1  # 奇偶性映射为 ±1
                            total_sum += parity * v          # 加权累加
                            sum += v
                    #print("当 cr1 全为0时，cr2 的计数统计:", cr2_counts)
                    #print("加权总和:", total_sum)
                    result_a = float(total_sum)/float(shots)
                    
                if all(bit == 'i' or bit == 'I' for bit in Ob[j]):
                    result_b = 1.
                    simulator = AerSimulator()
                    # Run and get counts
                    circ_b.measure(circ_b.qregs[0][0],  circ_b.cregs[0][0])
                    results = simulator.run(circ_b,shots=shots).result()
                    counts = results.get_counts(circ_b)
                    sum = 0
                    for k,v in counts.items():
                        cr1_part = k[:max_len]          # 前 m 位为 cr1
                        cr2_part = k[-1:]       # 后 n 位为 cr2
                        # 检查 cr1 是否全为0
                        if all(bit == '0' for bit in cr1_part):
                            sum += v
                    result_b = float(sum)/float(shots)
                else:
                    mat_rhob = np.zeros((1<<circ_b.num_qubits, 1<<circ_b.num_qubits),dtype= complex)
                    mat_rhob[0,0] = 1
                    circ_b_exe = circ_b
                    ob_str = list(Ob[j])
                    c2_num = 0
                    tag = 0
                    for o in range(len(ob_str)):
                        if ob_str[len(ob_str)-1-o] == 'x' or ob_str[len(ob_str)-1-o] == 'X':
                            circ_b_exe.h(o)
                            ob_str[len(ob_str)-1-o] = 'Z'
                        if ob_str[len(ob_str)-1-o] == 'y' or ob_str[len(ob_str)-1-o] == 'Y':
                            circ_b_exe.s(o)
                            circ_b_exe.h(o)
                            circ_b_exe.x(o)
                            ob_str[len(ob_str)-1-o] = 'Z'
                        if ob_str[len(ob_str)-1-o] == 'i' or ob_str[len(ob_str)-1-o] == 'I':
                            continue

                        circ_b_exe.measure(circ_b_exe.qregs[0][o],  circ_b_exe.cregs[0][tag])
                        tag += 1
                        c2_num += 1
                    #rhob = qiskit.quantum_info.DensityMatrix(mat_rhob)
                    #circ_b.h(0)
                    #rhob = rhob.evolve(circ_b_exe)
                    #print(ob_str,circ_b_exe)
                    #opb = SparsePauliOp.from_list([(''.join(ob_str), flag)])
                    #result_b = rhob.expectation_value(opb).real
                    simulator = AerSimulator()
                    # Run and get counts
                    results = simulator.run(circ_b_exe,shots=shots).result()
                    counts = results.get_counts(circ_b_exe)
                    c1_num = max_len
                    total_sum = 0
                    sum = 0
                    cr2_counts = {}
                    #print(counts)
                    for k,v in counts.items():
                        cr1_part = k[:c1_num]          # 前 m 位为 cr1
                        cr2_part = k[-c2_num:]       # 后 n 位为 cr2
                        # 检查 cr1 是否全为0
                        if all(bit == '0' for bit in cr1_part):
                            num_ones = cr2_part.count('1')
                            cr2_counts[cr2_part] = cr2_counts.get(cr2_part, 0) + v
                            parity = 1 if (num_ones % 2 == 0) else -1  # 奇偶性映射为 ±1
                            total_sum += parity * v          # 加权累加
                            sum += v
                    #print("当 cr1 全为0时，cr2 的计数统计:", cr2_counts)
                    result_b = float(total_sum)/float(shots)
                    #print("加权总和b:", result_b)
                result_i += result_a*result_b*o_coef[j]
                #print(o_coef[j],result_a,result_b)
            result += coef[i]*result_i
            #print(i,result)
        return  result


import random
import time
from qiskit import QuantumCircuit as QiskitQuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit_aer import AerSimulator
import numpy as np

def main():
    # 不同的采样次数设置
    shots_list = [1000, 2000, 3000]
    
    # 打开结果文件
    with open("knittingres.txt", 'w') as f:
        f.write("Qubits,Depth,Shots,Cutting_Time,Standard_Time,Cutting_Result,Standard_Result,Error\n")
        
        # 测试不同比特数（2到8，为了计算效率减少范围）
        for n in range(8, 9):
            # 测试不同深度（1到3）
            for d in range(1, 4):
                for shots in shots_list:
                    print(f"Testing n={n}, d={d}, shots={shots}")
                    
                    # 创建量子电路
                    qc = QuantumCircuit(n)
                    
                    # 构建电路
                    for i in range(n):
                        qc.add_gate(RX(i, random.uniform(0, 2*np.pi)))
                    
                    for j in range(d):
                        # 添加纠缠层
                        for i in range(n-1):
                            qc.add_gate(CZ(i, i+1))
                        # 添加单量子比特门层
                        for i in range(n):
                            gate_type = random.choice(['RX', 'RY', 'RZ'])
                            if gate_type == 'RX':
                                qc.add_gate(RX(i, random.uniform(0, 2*np.pi)))
                            elif gate_type == 'RY':
                                qc.add_gate(RY(i, random.uniform(0, 2*np.pi)))
                            else:
                                qc.add_gate(RZ(i, random.uniform(0, 2*np.pi)))
                    
                    # 创建随机可观测量
                    pauli_terms = []
                    coef_terms = []
                    # 添加单量子比特项
                    for i in range(n):
                        pauli_terms.append(f'Z {i}')
                        coef_terms.append(random.uniform(-1, 1))
                    # 添加两量子比特项
                    if n >= 2:
                        for i in range(min(3, n-1)):  # 最多添加3个两量子比特项
                            pauli_terms.append(f'Z {i} Z {i+1}')
                            coef_terms.append(random.uniform(-0.5, 0.5))
                    
                    O = Observable(pauli_terms, coef_terms)
                    
                    # 设置切割标签（将电路分为两部分）
                    split_point = n // 2
                    label = 'a' * split_point + 'b' * (n - split_point)
                    
                    # 使用电路切割方法
                    lce = largeCircuitExecution(qc, O, label)
                    
                    try:
                        # 执行电路切割
                        a, b, c = lce.CircuitCutting(label)
                        lce.ObservableCutting(label)
                        
                        # 计时电路切割方法
                        time_start = time.time()
                        cutting_result = lce.Execute(shots=shots)
                        cutting_time = time.time() - time_start
                        
                        # 使用qiskit标准采样方法计算
                        # 将自定义电路转换为qiskit电路
                        qiskit_circuit = custom_to_qiskit(qc)
                        
                        # 将可观测量转换为qiskit格式
                        qiskit_observable = custom_observable_to_qiskit(O, n)
                        
                        time_start_std = time.time()
                        standard_result = qiskit_sampling_expectation(qiskit_circuit, qiskit_observable, shots)
                        standard_time = time.time() - time_start_std
                        
                        # 计算误差
                        error = abs(cutting_result - standard_result)
                        
                        # 写入结果
                        f.write(f"{n},{d},{shots},{cutting_time:.4f},{standard_time:.4f},{cutting_result:.6f},{standard_result:.6f},{error:.6f}\n")
                        f.flush()
                        
                        print(f"  Cutting result: {cutting_result:.6f}")
                        print(f"  Standard result: {standard_result:.6f}")
                        print(f"  Error: {error:.6f}")
                        print(f"  Cutting time: {cutting_time:.4f}s")
                        print(f"  Standard time: {standard_time:.4f}s")
                        
                    except Exception as e:
                        print(f"  Error for n={n}, d={d}, shots={shots}: {str(e)}")
                        f.write(f"{n},{d},{shots},ERROR,ERROR,ERROR,ERROR,ERROR\n")
                        f.flush()

def custom_to_qiskit(qc):
    """将自定义电路转换为qiskit电路"""
    n_qubits = qc.qubit_number
    qiskit_circuit = QiskitQuantumCircuit(n_qubits)
    
    for i in range(1, len(qc.gate_list)):
        gate = qc.gate_list[i]
        qubits = gate.get_qubits()
        
        if gate.name == 'X':
            qiskit_circuit.x(qubits[0])
        elif gate.name == 'Y':
            qiskit_circuit.y(qubits[0])
        elif gate.name == 'Z':
            qiskit_circuit.z(qubits[0])
        elif gate.name == 'H':
            qiskit_circuit.h(qubits[0])
        elif gate.name == 'CX':
            qiskit_circuit.cx(qubits[0], qubits[1])
        elif gate.name == 'CZ':
            qiskit_circuit.cz(qubits[0], qubits[1])
        elif gate.name == 'RX':
            qiskit_circuit.rx(gate.get_para(), qubits[0])
        elif gate.name == 'RY':
            qiskit_circuit.ry(gate.get_para(), qubits[0])
        elif gate.name == 'RZ':
            qiskit_circuit.rz(gate.get_para(), qubits[0])
    
    return qiskit_circuit

def custom_observable_to_qiskit(observable, n_qubits):
    """将自定义可观测量转换为qiskit的SparsePauliOp格式"""
    pauli_list = []
    coeffs = []
    
    for pauli_str, coeff in observable.observable.items():
        # 创建Pauli字符串数组，初始化为'I'
        pauli_array = ['I'] * n_qubits
        
        # 解析Pauli字符串，如 'Z 0 Z 1'
        terms = pauli_str.split()
        for i in range(0, len(terms), 2):
            pauli_char = terms[i].upper()
            qubit_idx = int(terms[i+1])
            pauli_array[qubit_idx] = pauli_char
        
        # 转换为qiskit的Pauli字符串格式
        pauli_str_qiskit = ''.join(pauli_array)
        pauli_list.append(pauli_str_qiskit)
        coeffs.append(coeff)
    
    return SparsePauliOp(pauli_list, coeffs)

def qiskit_sampling_expectation(circuit, observable, shots):
    """使用qiskit采样计算期望值"""
    simulator = AerSimulator()
    
    # 对每个Pauli项分别测量
    total_expectation = 0.0
    
    for pauli, coeff in zip(observable.paulis, observable.coeffs):
        # 复制电路
        meas_circuit = circuit.copy()
        
        # 添加基变换以测量该Pauli项
        for i, pauli_char in enumerate(pauli):
            if pauli_char == 'X':
                meas_circuit.h(i)
            elif pauli_char == 'Y':
                meas_circuit.sdg(i)
                meas_circuit.h(i)
        
        # 添加测量
        meas_circuit.measure_all()
        
        # 运行电路
        result = simulator.run(meas_circuit, shots=shots).result()
        counts = result.get_counts()
        
        # 计算该Pauli项的期望值
        pauli_expectation = 0.0
        for bitstring, count in counts.items():
            # 计算该比特串对应的特征值（+1或-1）
            eigenvalue = 1
            for i, pauli_char in enumerate(pauli):
                if pauli_char != 'I':
                    # qiskit的比特顺序是反的
                    bit_index = circuit.num_qubits - 1 - i
                    if bitstring[bit_index] == '1':
                        eigenvalue *= -1
            pauli_expectation += eigenvalue * count
        
        pauli_expectation /= shots
        total_expectation += coeff * pauli_expectation
    
    return total_expectation

if __name__ == "__main__":
    main()